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ABSTRACT 

This report is subdivided in three parts. The first one reviews a new approach to the computation of 
inviscid flows using Cartesian grid methods. The crux of the method is the curvature-corrected symmetry 
technique (CCST) developed by the present authors for body-fitted grids. The method introduces ghost 
cells near the boundaries whose values are developed from an assumed flow-field model in vicinity of the 
wall consisting of a vortex flow, which satisfies the normal momentum equation and the non-penetration 
condition. The CCST boundary condition was shown to be substantially more accurate than traditional 
boundary condition approaches. This improved boundary condition is adapted to a Cartesian mesh 
formulation, which we call the Ghost Body-Cell Method (GBCM). In this approach, all cell centers 
exterior to the body are computed with fluxes at the four surrounding cell edges. There is no need for 
special treatment corresponding to cut cells which complicate other Cartesian mesh methods. The merits 
of the procedure are indicated by the computation of the compressible flow about basic circular cylinders 
as well as of airfoil applications and detailed comparisons to body-fitted grid computations. The results 
show the surface non-penetration condition to be satisfied in the limit of vanishing cell size and the 
method to be second-order accurate in space. In addition, the results are independent of the position of 
the body with respect to the grid. The flexibility of the approach is illustrated with an application to 
the flow over a multiple-body configuration. 

The original version of the Ghost Body-Cell Method, had the typical Cartesian grid requirement that 
any grid clustering near the body must be maintained to the far field boundary. This is not a significant 
penalty in 2D calculations, but it may become significant in 3D calculations. As a consequence, in the 
second part of the paper, we have introduced a far-field grid coarsening, based on the iblanking approach, 
which maintains the structured nature of the grid while computing only the active cell centers. The 
far-field grid coarsening has been applied to subsonic and transonic flows about an NACA 0012 airfoil, 
indicating that the results are practically unaffected by the coarsening. Meanwhile, the number of the 

computed cells is reduced to 40% 50% and the CPU time to 34% 47% of the original CPU time, 

depending on the particular test case. 

The GBCM on far-field coarsened Cartesian grids appears to be particularly well-suited for design 
optimization problems. The third part of the paper is dedicated to this issue and presents an efficient for- 
mulation for the robust inverse design of compressible fluid flow problems on far-field coarsened Cartesian 
grids. The design approach has five essential features: First, the procedure utilizes the GBCM on coars- 
ened Cartesian grids. Second, approximate but efficient design sensitivities are obtained using a discrete 
adjoint formulation, based on an artificially dissipative flow solver in order to obtain robust sensitivity 
derivatives in the presence of potentially noisy or non-smooth objective functions. The third feature is the 
implementation of quasi-time-dependent adjoint equations. Next, the procedure involves what we term 
progressive optimization, whereby a sequence of operations is performed containing a partially converged 
flow solution, followed by a partially converged adjoint solution followed by an optimization step. The 
final feature is the use of progressively finer grids for the progressive solution of the flow field, while only 
using a reduced domain for the solution of the adjoint equations. 

The inverse design of airfoils in transonic flow conditions indicated that the method is accurate, 
robust and highly efficient, with converged design optimizations produced in no more than the amount 
of computational work to perform two flow analysis on the finest mesh. In addition, a test case involving 
inverse design of a NACA 0012 airfoil in subsonic flow conditions showed that the approach is able to 
bring the convergence of the optimization process to machine zero. 

‘Professor, Dipartimento di Ingegneria Meccanica e Gestionale, Sezione Macchine ed Energetica. 

TVice President, Education & Outreach. Also Prof. Aero, h Ocean Eng., Virginia Tech. 
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INTRODUCTION 


The continuous growing power of computers is encouraging engineers to rely on computational fluid 
dynamics (CFD) for the design and testing of new technological solutions. Numerical simulations allow 
the analysis of the phenomena without resorting to expensive prototypes and complex experimental mea- 
surements. On the other hand, while simple geometries discretized by regular grids are efficiently handled 
by codes and hardware, flows with complex geometries requiring body-fitted grids are still challenging 
problems for today computers. The grid generation process, whether for composite block structured 
grids or unstructured grids is often laborious and is not always as automated as desired. The conceptual 
design process often requires detailed body-fitted grids over wide classes of different geometric topologies. 
Grid generation for complex topologies often takes weeks or even longer to develop and this becomes a 
major impediment to using CFD in design. Further complications may occur due to moving boundaries 
or problems which require continuing regeneration or deformation of grids. 

One approach to develop a simplified alternative to complex grid generation involves the use of 
Cartesian grid methods. In these methods the body does not coincide with grid surfaces or grid volume 
edges. Instead, the body is immersed in the grid and the complication of generating a body-fitted grid 
is exchanged for complications associated with the implementation of the surface boundary conditions. 
Furthermore, since the body surfaces are not aligned with the grid it is nearly impossible to adequately 
resolve high Reynolds number boundary layers with Cartesian mesh methods. Nonetheless, there are 
many classes of flow problems which benefit from the efficiency of Cartesian mesh approaches. 

Cartesian mesh methods have their origin in transonic potential flow computations 1,2 . Early Cartesian 
Euler computations were performed in Refs. 3-5. More recent computations 6-8 involved the use of 
grid adaptation and have produced very good results on some very complex configurations. Aftosmis 
and Berger 9 have also shown that Cartesian mesh computations may be characterized by a truncation 
error which is approximately one order of magnitude smaller than truncation error of computations on 
triangular or quadrilateral body-fitted grids. 

For completeness, we mention two other approaches which hold promise for avoiding grid generation 
difficulties. One is a Cartesian grid method which simulates viscous flows using vorticity confinement 10 
and the other is an approach in which a body force term is used to simulate moving boundaries 11,12 . 

Cartesian grid finite-difference schemes for computing fluid dynamic problems, e.g ., Ref. 5, have 
proven to be quite efficient, but require a complex topological discrimination of the cells located near the 
body profile. Cartesian grid finite-volume methods, e.g. Ref. 13, are more straightforward in that they 
allow to use composite elements (cut cells) determined by the intersection between the Cartesian grid 
and the body profile. On the other hand, these methodologies may suffer stability problems when such 
composite elements become very small. 

The use of Cartesian grids without composite elements avoids this stability problem. On the other 
hand, it significantly increases the problem of enforcing the non-penetration condition at solid walls 
immersed in the Cartesian grid. For inviscid flows, the simple non-penetration condition can be easily 
enforced at solid walls if body-fitted grids are used, because the grids are aligned with the wall. For 
Cartesian grids, appropriate flow conditions must be enforced at the cell centers surrounding the body. 
The accuracy of the computed flow solution obviously depends on how accurately the flow conditions 
are enforced at these cell centers. As a consequence, an accurate model of the flow near the solid wall is 
required if accurate flow solutions on Cartesian grids are sought. 

Recently, a model able to accurately represent inviscid flows near the solid walls has been introduced 
by the present authors 14 . It has been named CCST ( curvature- corrected symmetry technique ) and has 
been developed for body-fitted grids. The method considers solid walls as boundaries immersed in the flow 
field and enforces conditions at ghost cell centers located inside the body in a position close to the wall. 
The ghost cell values are developed from an assumed flow-field model in the vicinity of the wall consisting 
of a vortex flow, which satisfies the normal momentum equation and the non-penetration condition. This 
flow-field model enforces symmetry conditions for entropy and total enthalpy along a normal to the body 
surface. Then, the flow computations are performed using these ghost cell centers without taking into 
account the presence of the wall. Obviously, the fluid dynamic parameters are continuously updated 
at these ghost cell centers in order to enforce the wall impermeability condition at each stage of the 
computation. 

The CCST boundary condition has been shown to be significantly more accurate for a variety of 
inviscid flows and has been shown to be at least one order of magnitude more accurate in inviscid two- 
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dimensional flows 14 as well as in inviscid three-dimensional flows 15 . Other researchers have successfully 
implemented this boundary condition, e.g. Ref. 16 have reported excellent results using the COST in a 
preconditioned Euler code. 

The first part of the present paper is dedicated to the review of an innovative and accurate numerical 
technique to enforce non-penetration conditions at solid walls in Cartesian grid computations. This 
new approach utilizes the curvature corrected symmetry technique applied to Cartesian grids, which we 
call the Ghost Body-Cell Method (GBCM). This method, in some of our earlier papers, was called the 
Immersed Body Methodology 17,18 . This approach uses only rectangular cells without having to define cut 
cells in the vicinity of the solid boundary. We first describe the governing equations and the classical 
boundary conditions to enforce the non-penetration condition at solid walls for body-fitted grids. Then the 
improved boundary condition called the curvature-corrected symmetry technique (CCST) is presented. 
Finally, the CCST is extended to model the flow near solid walls immersed in a Cartesian grid, GBCM. 
Details related to sharp corners and trailing edges are also described. The merits of the procedure are 
indicated by the computation of the compressible flow about basic circular cylinders as well as of airfoil 
applications and comparisons to body-fitted grid computations. The flexibility of the method is illustrated 
with an application to the flow over a multiple-body configuration. 

In the aforementioned applications, the method is shown to produce accurate results for elementary 
two-dimensional geometries with relatively coarse grids. Unfortunately, in usual Cartesian grid formu- 
lations, any grid clustering near the body must be maintained to the far field boundary. This is not a 
significant penalty in 2D calculations, but it becomes significant in 3D calculations. As a consequence, 
the second part of the paper is dedicated to far-field grid coarsening to address this issue and improve 
the method’s efficiency. Ref. 19 provides a basic approach to local refinement of structured grids. The 
suggested procedure is based on the idea of iblanking. The iblanking concept was introduced in Ref. 20 
as a device to insert geometry into a structured grid by extending the grid inside the body and then 
blanking out the interior portion. In Ref. 19, the iblanking algorithm is extended to a coarsened struc- 
tured grid, so that the blanked region is skipped, requiring neither storage of variables nor solution of 
equations. In such a way, incompressible Navier-Stokes equations in Ref. 19, were efficiently solved using 
a finite-difference method on a Cartesian grid. In the second part of the present paper, the iblanking 
approach is extended to the present finite-volume formulation. Results for airfoils will be presented with 
and without far-field coarsening and compared in terms of accuracy and efficiency. 

The Ghost Body-Cell Method on coarsened Cartesian grids appears to be particularly well-suited 
for design optimization problems. The third part of the paper is dedicated to the presentation of the 
efficient design procedure that we utilize. This approach was originally developed for body-fitted grids 21 
and is adapted here for Cartesian grids. It has five essential ingredients. First, the procedure utilizes 
an accurate finite-volume method on coarsened Cartesian grids (GBCM) in order to accurately evaluate 
the objective function. Next, approximate design sensitivities are efficiently obtained using a discrete 
adjoint formulation, based on a simplified artificially dissipative flow solver in order to obtain robust 
sensitivity derivatives in the presence of noisy or non-smooth objective functions. This adjoint solver 
is only used for the solution of robust sensitivities. The third ingredient of the procedure is the use of 
quasi-time-dependent adjoint equations. The fourth ingredient of the procedure involves a progressive 
optimization, whereby a sequence of operations, containing a partially converged flow solution, followed 
by a partially converged adjoint solution followed by an optimization step is performed. Both the flow 
and the adjoint equation solutions are performed on Cartesian grids with body wall boundaries immersed 
in the Cartesian grid. Furthermore, the adjoint equations are solved on a reduced domain compared to 
the one used to compute the flow solution. The fifth ingredient is the use of progressively finer grids 
for the progressive solution of the flow field, while only using a reduced domain for the solution of the 
adjoint equations. In the third part of the paper the optimization procedure will be presented together 
with some computed results for the inverse design of airfoils. 

1. GHOST BODY-CELL METHOD 


1.1 Governing Equations 

The vector form of the two-dimensional Euler equations in a Cartesian coordinate system can be 
written as: 
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0 , 


( 1 ) 


q f + f r + g y = 

where the vector of conserved variables is given by: 

q = [p, pu, pv, pe 0 ] T , (2) 

and the inviscid flux vectors are: 


f = [pu, p + pu 2 , puv, puh 0 ] T , 

g = [pv, puv, p + pv 2 , pvh 0 ] T , 


( 3 ) 


where p, p, u, v are the pressure, the density and the two Cartesian components of the velocity vector, 
respectively, while eo and ho represent the total energy and the total enthalpy per unit mass, respectively. 

In the work presented here, a quadrilateral, cell-centered finite- volume method has been implemented 
with a flux-difference-splitting algorithm 22 as the flow solver. A MUSCL-type extrapolation, with a 
formal second-order accuracy in space, has been applied to extrapolate the physical variables onto the 
left and right sides of each cell edge. However, the methods developed in this paper is applicable to other 
flow solvers. 


1.2 Wall Boundary Conditions 

Classical non-penetration condition 

We consider a cell-centered finite- volume representation of Eqs. (1) and focus our attention on a 
quadrilateral cell of the computational volume. The flux vector across a cell edge can be expressed as: 

F = [pv, n x p + puv, n y p + pvv, ph 0 v] T , (4) 

where n x and n y are the direction cosines of the normal to the cell edge and v is the corresponding 
velocity component: 


v = n x u + n y v. (5) 

Assume that this cell face is located on a solid boundary so that v = 0. We see from Eq. (4) that 
the only quantity required to compute the flux vector at this cell edge is the pressure at the wall. For 
a cell-centered scheme, there are several conventional ways to evaluate this pressure. The most common 
approach is to extrapolate from the interior of the flow field. A first-order pressure extrapolation is often 
used whereby the wall pressure is taken as the value associated with the nearest cell center. A second- 
order pressure extrapolation (P-II) is also sometimes used with the wall pressure linearly extrapolated 
from the two nearest cell centers: 


Pw = -jpi ~ 2^2, (6) 

where subscripts 1 and 2 refer to these cell centers (Fig. 1). In another approach, the pressure may 
be extrapolated to satisfy the normal momentum equation. Yet another boundary condition procedure 
utilizes a characteristic boundary condition. A description and comparison of these techniques together 
with their application to a shock reflection problem are reported in Ref. 14. 


CCST: a more accurate non-penetration condition 

The curvature corrected symmetry technique (CCST) described in Ref. 14 utilizes ghost points and 
a vortex model in the vicinity of the body surface. The procedure consists in first locating two rows 
of cell centers closest to the body surface inside the flow domain (points 1 and 2 in Fig. 1). Then two 
extra rows of cells are inserted (image cells or ghost cells) outside the computational flow field which are 
created by reflecting the internal cells 1 and 2 symmetrically with respect to the body surface (points — 1 
and —2 in Fig. 1). The fluid dynamic variables at the ghost cell centers are utilized to implement the 
non-penetration condition at the wall. The ghost cell values are developed from an assumed flow-field 
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model in vicinity of the wall, consisting of a vortex flow with locally symmetric distribution of entropy 
and total enthalpy along a surface normal. This flow model satisfies the normal momentum equation, the 
non-penetration boundary condition and enforces antisymmetric normal derivatives dS/dn and dho/dn 
along the normal to the body surface in the vicinity of the body. These entropy and total enthalpy 
distributions will produce zero normal derivatives when the flow is irrotational and will satisfy Crocco’s 
relationship even when vorticity is present near the surface. The method has shown to produce superior 
accuracy 14 compared to the classical surface boundary condition procedures discussed in the previous 
subsection. 

The pressure at the image cell centers is obtained from the integration of the normal momentum 
equation (for a steady inviscid flow): 


11 ■ 

P-i = Pi - Pi-ET^nt, (7) 

ilj 

where the subscript i is equal to 1 or 2 , An,; indicates the distance between the cell centers (*) and (— i), 
and Uj is the velocity component tangential to the body: 


u = n v u—n x v. ( 8 ) 

Note that 7?,; is the signed local radius of curvature of the wall, positive if the center of curvature is on 
the same side of the body, negative if the center of curvature is on the side of the flow field. 

The density and the velocity component tangential to the body are obtained from our vortex model: 


P-i 

= Pi 

V-i 

= u] + 


P-i 

Pi 


1/7 


0 ) 


27 

7-1 


El 

Pi 


p-i 

p-i 


. -2 -2 
+ '•, -V-i, 


where the subscript i is equal to 1 or 2 and 7 is the constant ratio of specific heats. 

The velocity component normal to the body is obtained from the non-penetration condition, following 
Ref. 14: 


V-l 

v - 2 


-Vl , 

3fj_i + 

2 [fq + (Ui 


{io)/2] 1 


I Pi + (pi - P2)/2 
P-1 + (p-i - P-2 )/2 


( 10 ) 


Numerical experiments have also shown that V -2 can be approximated by a value antisymmetric to ib • 
Then Eqs. (10) can be simplified as follows: 


v-i = ~Vi, ( 11 ) 

where the subscript i is equal to 1 or 2 . 

In Ref. 14 we show dramatic advantages of the CCST over conventional methods with regard to 
numerical entropy generation, total pressure loss and grid convergence. 

1.3 Cartesian Grid Method 

In this section, the curvature corrected symmetry technique (CCST) will be extended to non-body 
fitted Cartesian grids. The method is developed for two-dimensional flows, but may be extended to three 
dimensions. The new numerical technique, the Ghost Body-Cell Method, is described by the following 
steps: 

1. The body is defined by a dense, ordered sequence of ( x,y ) coordinate pairs. The definition of the 
surface is completely independent of the mesh. 

2. A Cartesian grid is selected and the computational flow field is subdivided into finite-volume cells. 
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3. The cell centers are discriminated between those located at the interior of the body and those 
located externally. The internal cell centers are flagged with A = 0, while the external ones are 
flagged with A = 1. 

4. The row of cell centers closest to the body and located at the interior of the body is singled out. 
These cell centers are indicated by circles in Fig. 2 which presents a part of the cell center net and 
the curved body profile (the right-hand side of the solid line is interior to the body). 

5. A second row of cell centers located at the interior of the body and adjacent to the first row is 
singled out. These cell centers are indicated by triangles in Fig. 2. 

6. For each of the interior cell centers of the two rows described above (ghost cell centers), a correspond- 
ing symmetric point is determined at a location exterior to the body and reflected symmetrically 
with respect to the body surface. To locate these symmetric points, we first determine the straight 
line passing through the cell center and orthogonal to the body. On this line, the point symmetric 
of the interior cell center with respect to the body is determined. These symmetric points are 
indicated by dark squares in Fig. 2. As an example, the point B and the interior cell center A in 
Fig. 2 are symmetric with respect to the body. 

7. The four cell centers surrounding each symmetric point are determined. As an example, points C, 
D, E, F are the cell centers surrounding point B in Fig. 2. Note that one or even two of these four 
cell centers may be located inside the body. 

8. The value of the conserved variables q at the point B are determined by a bilinear interpolation: 


qs = ki(x - x c )(y ~ yc) + h(x - x c )+ n9 x 

h(y-yc) + k 4 (U) 

where the constants kj (i = 1, . . .4) are determined by the values of the conserved variables q at 
the cell centers C, D, E, F. The bilinear interpolation is also used in situations when one or two of 
the cell centers surrounding point B are located inside the body surface. The internal values are 
determined from the CCST boundary condition (described next). 

9. The values of the physical variables at the cell centers belonging to the previously defined two rows of 
internal cell centers are determined according to the CCST boundary condition, Eqs. (7), (9), (11). 
As an example, the physical variables at the cell center A are given by: 


PA = Pb — Pb 


if*" 


(13) 


Pa = Pb 



(14) 


-2 _ -2 . 2 7 (PB PA\ 


(15) 


VA = ~VB (16) 

In Eqs. (13)— (16), An is the distance between points A and B, while R is the signed local radius of 
curvature of the body. 

10. The time-dependent computation of all conserved variables at cell-centers located external to the 
body now proceeds without any other boundary condition. Each cell center utilizes the fluxes 
at surrounding cell edges without further approximation. Cell centers located near the boundary 
utilize the interior (ghost) cells. 
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The time dependent computation of the flow field starts from the undisturbed flow condition enforced 
at all the cell centers located outside of the body, while stagnation flow conditions are enforced inside 
the body. In addition, the computed residuals are multiplied by A. These assumptions are irrelevant 
to the flow computation and have been used in order to avoid any computational problem at the cell 
centers located inside the body. In addition, the last assumption enables a check of the variations of the 
residuals, which would never converge if the residual of the cell centers located inside the body are taken 
into account. 

A special procedure based on double-valued cell centers is required near unresolved thin surfaces and 
sharp corners, such as at the sharp trailing edge of an airfoil. At sharp corners, one cell center inside the 
geometry may be the ghost cell center for one side of the corner surface as well as for the other side. Also, 
a ghost cell center pertaining to one side of the corner surface may be located in the flow field on the 
other side of the sharp corner. One typical situation is shown in Fig. 3, where the first row of cell centers 
are again represented by circles and the second row by triangles. Note that point R has a symmetric 
point on the upper side of the trailing edge and a second symmetric point on the lower side. A similar 
situation holds for point U. In addition, point S belongs to the second row of ghost cell centers and has a 
symmetric point on the lower side of the trailing edge, while cell center T has a symmetric point located 
on the upper side. As a consequence, the cell centers R, S, T, U will be double-valued. The procedure 
we use in this region is described by the following steps: 

1. The sharp corner region is defined by detecting the body region where at least one of the ghost 
cell centers pertaining to one side of the sharp corner does not differ from the ghost cell centers 
pertaining to the other side or is located inside the flow field. In Fig. 3, this region includes points 
R, S, T, U. Note that cell centers V, W, Q, Z do not belong to this region because the first two 
have symmetric points on the lower side and the other two cell centers have symmetric points on 
the upper side so that no cell center is double- valued. 

2. The fluid dynamic properties at the ghost cell centers are computed by applying the previous 
procedure. The cell centers R and U, located inside the body, are characterized by two different 
sets of fluid dynamic properties evaluated by applying the previous procedure to both sides of the 
trailing edge. Also, the cell centers S and T, located in the flow field, are characterized by a first set 
of properties corresponding to the flow field and by a second set computed by applying the previous 
procedure to the lower and to the upper side of the trailing edge, respectively. 

3. After computing the fluxes in a traditional manner, we recompute the fluxes involving the double 
values at cell centers R, S, T, U. 

4. We recompute the fluxes downstream of the trailing edge, involving the points R and U, by averaging 
the two sets of values at these points. 

In GBCM, all cell centers exterior to the body are computed with fluxes at the four surrounding cell 
edges. There is no need for special treatment corresponding to cut cells which complicate other Cartesian 
mesh methods. Special treatment for the evaluation of ’’ghost cell centers” as well as surface normals 
and radii of curvature are required. Since the body is represented pointwise and is independent of the 
mesh, it may be easily defined by a sufficient point density to allow straightforward accurate numerical 
calculation of the surface normals and radii of curvature. For the present test cases, the description of 
the body profile with 200 points, connected by segments, produces a sufficiently accurate calculation of 
An and R. 

1.4 GBCM Results 

First, the subsonic compressible inviscid flow over a circular cylinder has been analyzed to evaluate 
the merits of the Ghost Body-Cell Method for Cartesian grids. A free-strearn Mach number of = 0.38 
has been prescribed. This flow condition produces a maximum Mach number of approximately 0.92 on 
the surface of the cylinder. The computations have been performed using three Cartesian grids: 110 x 110, 
56 x 56, and 28 x 28. The numbers of cell centers located in the first row surrounding the body on the 
internal side were 128, 64 and 32, respectively. The far-field boundary was located at approximately 20 
diameters. 
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Three parameters have been used to evaluate the effectiveness of the Cartesian grid technique. These 
include the Lo-norm of the normalized surface velocity component normal to the body, {v &), the Mach 
number at the leading edge, M) e , and the Mach number at the trailing edge, M te . These three parameters 
are plotted in Fig. 4 versus 1/N' 2 , where N is the number of cell centers surrounding the body on the 
internal side. In order to generate the data for the curves in Fig. 4, the three parameters have been 
computed by interpolating the conserved variables at the body points from the four surrounding cell 
centers. Note that one or even two of the four cell centers may be located inside the body. Fig. 4 
indicates that GBCM enforces the non-penetration condition in the limit of vanishing cell size, and 
appears to be second-order accurate in space. The errors in Fig. 4 are very small (of the order of 10 -4 ) 
when the finest grid is employed. Note that the free-stream velocities are scaled to 1. 

The flow has been also computed using a body-fitted polar grids with 128 x 32 cells, with the larger 
number corresponding to the number of cells located in the circumferential direction. The wall non- 
penetration condition has been enforced with a classical second-order pressure extrapolation technique 
(P-II). The number of cell centers surrounding the body on the internal side is the same for polar grid 
and the corresponding Cartesian grid. The results with the Cartesian grids (using GBCM) and with 
the body-fitted polar grids have been compared. Figs. 5 and 6 present the contours of the normalized 
total pressure. The contours of total pressure found using CCST methodology with the Cartesian-grid 
technique are very regular, with a very small total pressure increase, corresponding to a maximum error 
equal to 0.2%, quite close to the maximum error found with the CCST on polar grids (see Ref. 14). 
In contrast, the P-II technique yields significant total pressure errors, with a maximum equal to 1.4%. 
Figs. 5 and 6, together with other results reported in Ref. 17, indicate that the use of the Ghost Body- 
Cell Method on Cartesian grids preserves the superior accuracy of the curvature-corrected symmetry 
technique on body-fitted polar grids and shows dramatic advantages with respect to the widely used 
second-order pressure extrapolation technique on body-fitted polar grids. 

Next, more general two-dimensional airfoil geometries have been considered. Airfoils are of course 
an important class of geometries, but they also serve as a test of the method for sharp corners, such 
as near sharp trailing edges. For these cases, the double-valued ghost-cell-center procedure for sharp 
corners, outlined at the end of §1.3, has been applied. Results have been obtained for a NACA 0012 
airfoil. Various Mach numbers and angles of attack, corresponding to subsonic as well as transonic flow 
conditions have been considered. Here we present one of the transonic results, corresponding to a free- 
stream Mach number l/ x , = 0.8 at an angle of attack equal to 1.25°. Fig. 7 presents the computed Mach 
number distribution on the airfoil surface (continuous lines) compared with the body-fitted grid Mach 
number distribution (symbols) obtained using the CCST technique. The two sets of results in Fig. 7 
agree quite well. The computed lift and drag coefficients are: Cl = 0.3535, Cd = 0.0232. These values 
compare quite well with the results obtained by Salas and South (as discussed in Ref. 23, where their 
values Cl = 0.3474, Cd = 0.0221 are described as accurate. 

The influence of the airfoil leading edge position with respect to the grid has also been investigated. 
A transonic test case with = 0.85 at zero incidence has been considered. The flow has been first 
computed locating the leading edge precisely at a cell center, which will be considered as the reference 
test case in the following. Then four more computations have been performed by moving the leading edge 
by —50%, —25%, 25%, 50% of the local cell dimension in the x direction. Fig. 8 shows the Mach number 
distribution in the leading-edge region for three of these computations: the continuous line corresponds to 
the leading edge located coincident with the cell center, the circles indicate the results corresponding to 
the leading edge located 50% upstream, and finally the X symbols correspond to the leading edge located 
50% downstream. Note that the abscissa is measured from the leading edge position. The position of the 
airfoil profile with respect to the cell center net is also shown in the small window in Fig. 8, using the same 
symbols already used for the Mach number distribution. Note the position change of many cell centers 
when changing the location of the leading edge with respect to the cell center net: many cell centers 
located inside the body are moved to the flow field and vice versa so that the geometrical configuration 
presents large changes. Nevertheless, the three sets of Mach number distributions coincide extremely well 
despite the effects of the curvature in this region and the changes in the geometrical configuration. This 
finding conclusively indicates that the method is practically independent of the airfoil profile location with 
respect to the grid. In order to reenforce this conclusion, Table 1 presents the mean square difference 
between the surface values corresponding to the different grid locations for some fluid-dynamic properties 
(pressure, density, Mach number, and the two velocity components). Also shown in Table 1 is the mean 
square error of the normalized total temperature on the body surface, (T° — 1), and of the normalized 



velocity component normal to the body surface, (v), together with the lift, Cl, and drag coefficients, 
Cd, for the five prescribed locations of the leading edge. Note that Table 1 contains interpolated values 
which are used only for assessment purposes. Table 1 quantifies the effects of moving the leading edge 
with respect to the grid and establishes that the results are practically independent of the airfoil profile 
location with respect to the grid. This table also indicates that the normalized total temperature errors 
are very small (of the order of 10 -3 ) as well as the normalized velocity components normal to the body 
(of the order of 10 -4 ). This result is important in order to effectively utilize this method for application 
involving design optimization. 

In order to show that the Ghost Body-Cell Method can be easily extended to multiple body configu- 
rations, the external/internal flow over the so called bi-NACA 0012 double airfoil is considered. The test 
case has been suggested for the GAMM workshop 24 on the numerical simulation of compressible Euler 
flows held at INRIA on June 1986. The geometry of this test case, plotted in Fig. 9, is made up of two 
NACA 0012 airfoils separated by a distance equal to 50% of the airfoil chord. Three flow conditions 
have been suggested in Ref. 24: a subsonic flow condition, a transonic flow with an internal shock, and 
a transonic flow with internal and external shocks. This last flow condition corresponds to an upstream 
Mach number = 0.55 at an angle of attack of 6° and has the most interesting flow features. Indeed, 
on the upper side of both airfoils, there is a supersonic expansion due to the curvature of the leading 
edge, terminated by a sudden recompression. Also, the internal region between the two airfoils is shaped 
as a nozzle where a supersonic expansion terminated by a shock takes place. Five contributors of Ref. 
24 presented results for this case in detail and one author provided only the corresponding lift and drag 
coefficients. 

The computation has been performed using a Cartesian grid with 140 x 202 cells with 104 cell centers 
located in the first row surrounding the body on the internal side. The corresponding constant Mach 
number contours are plotted in Fig. 9, while the computed lift and drag coefficients are reported in 
Fig. 10 where they are compared to the values computed by the six contributors 24 . Fig. 10 shows the lift 
coefficient versus the drag coefficient for the lower and the upper airfoils: the present results are indicated 
by dark symbols while the other results from Ref. 24 are indicated by light symbols. We see from Fig. 10 
that the present results are in general agreement with those of Ref. 24 and that they are quite close to 
the results by Couallier et al. (right triangles), which were considered to be the most accurate 24 . 

2. FAR-FIELD COARSENING 


2.1 Iblanking Approach 

The iblanking approach 19 is used here to preserve the structured nature of the Cartesian grid. All 
the cell centers are conserved, even in presence of grid coarsening, but only the active (non-iblanked) 
cell centers are computed. Consider Fig. 11, where the active cell centers are indicated by solid squares, 
while the iblanked cell centers are indicated by open squares. Before the grid coarsening, four cell centers 
located at j\, j\ + 1, j\ + 2 and j\ + 3 are active at i\ as well as at i i + 1. After the grid coarsening 
has been applied, four cell centers are still active at i\, while only two cell centers are active at i\ + 1, 
because those located at j\ + 1 and at j\ + 3 have been iblanked, i.e., they are still present in the data 
structure, but the equations are not integrated at these cell centers. Note that the still active cell centers 
at (ii + 1, j\) and (*i + 1, j\ + 2) may be considered to have been geometrically moved to the positions 
indicated by dark circles. Obviously, the equations will be evaluated at these active cell centers and will 
be performed with the correct geometric data. In these computations, data will be taken from active cell 
centers, i.e., the computation at (ii + 1, j\) will take data from (i i + 1, j\ -f 2), while skipping the data 
at (ii + 1, j\ + 1), because it is not active. 

The criterion for iblanking the cell centers is the following: First, the computational flow field is 
subdivided into five regions. The first rectangular region is located at the center of the computational 
flow field and includes the body (refer to Fig. 12); no coarsening is performed in this region in order 
to preserve an adequate dense mesh close to the body. The other four regions are located upstream 
and downstream of the body in the x and y directions, respectively as indicated in Fig. 12. Consider 
the regions upstream and downstream of the body in the x-direction (regions U and D in Fig. 12). As 
you move away from the body, Ax,; increases due to grid stretching. In these regions, far upstream and 
downstream of the body, it is clear that Ax, will be very much greater than A yj. A rational criterion 
for grid coarsening in these regions is that the grid will be coarsened when at a given x,; we find that 
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A yj + Ayj + i < A Xi. Similarly in the regions above and below the body (regions A and B in Fig. 12), the 
grid stretching has A yj increasing as you move away from the body. In these regions a rational criterion 
is to coarsen the grid when at a given yj we find that A Xi + Ax,;+i < Ayj. In each of the regions, the 
coarsening is applied sweeping from the center of the flow field to the corresponding outer boundary. 

At this point we illustrate the grid coarsening procedure. At the beginning, all the flow field cell 
centers are flagged with 1, then we start the grid coarsening in one of the four coarsening regions U, 
D, A, B of Fig. 12. Let us suppose to start the grid coarsening in the region downstream of the body 
(D) where the coarsening starts at the lowest i value. Referring to Fig. 11, we suppose that the grid 
coarsening begins with ij + 1. In Fig. 11, the cell centers after grid coarsening are indicated by solid 
circles. The computed flow quantities will be stored in the corresponding solid squares (below). The 
open squares will be treated by iblanking. The way that this is implemented is as follows: we sweep 
along the line i = i\ + 1 and check if two nearby cells coalesce, (indicated by j\ and j\ A 1 in Fig. 11) 
and determine a new cell characterized by a Ay <0.9 A;r. If this is true, the two cells coalesce to form 
a new larger cell at i\ -F 1 by merging the j\ and j\ + 1 cells. The merged cells are extended to the outer 
boundary for all i > i\ -F 1. Meanwhile all the cell centers j\ + 1 are iblanked for i > + 1 and flagged 

with a 0. Obviously the cell centers j\ for i > i\ + 1 remain active and flagged with 1, but with updated 
geometric properties. Sweeping the whole i i + 1 line complete the grid coarsening at i\ -F 1. Then we 
move to i\ + 2 and we repeat the previous process by considering only the cell centers flagged with 1, 
because those flagged with 0 have already been iblanked. A similar process of grid coarsening is applied 
up to the outer boundary and repeated in the other three regions, so that a far-field coarsened grid is 
finally obtained. 

Next, in order to solve the equations at each active cell center, we need to determine the nearby active 
cells. This may be simply performed by sweeping all the i or j lines and determining the distances, in 
terms of i and j, between the active cell center and the nearest active cell centers. In such a way, the values 
of Ai + ,Ai~ , A j + , A j~ can be evaluated. Obviously, in absence of any coarsening, Ai + = A j + = 1, 
A i~ = A j~ = —1, while the absolute values of these quantities may be larger in coarsened regions. Then, 
in order to solve the equations only at each active cell center, the previous sweep is repeated and the i 
and j coordinates of the active cell centers (flagged with 1) are inserted in two one-dimensional arrays 
and the computations are performed only at these coordinates. In addition, the cell centers located inside 
the body are flagged with a 0, so that these cell centers are not computed. 

Some particular treatment is required at the coarsening elements. In Fig. 11, a coarsening element is 
the zone including the cells (i± r j i), (ii,ji + 1), (i i -F 1, j\). Typically, a coarsening element contains two 
small cells and one adjacent larger cell, obtained from the coalescence of two smaller cells. The fluxes at 
the vertically-oriented cell edges cannot be computed in a straightforward way because of the iblanked 
cell center (i i + 1, ji -F 1) and because of the change of the geometric definition of the active cell center 
(ii + 1, j\). A typical topological picture of the region surrounding a coarsening element in the coarsening 
region downstream of the body is reported in Fig. 13. In this figure, the dark circles indicate cell centers, 
some of which are designated by a capital letter. The X symbols and the corresponding numbers indicate 
cell edges and the + symbols and the corresponding numbers locate points used for extrapolation of the 
left and right states. Note that the scales of A" and Y in Fig. 13 are different. As a consequence the cell 
dimensions are different from the representation in Fig. 11, e.g., the cell P is almost a square while in 
Fig. 13 it appears to be a rectangle. The left and right states at the cell edges labelled with the numbers 
1, 5, 9, 13, 14 can be evaluated by means of the traditional linear extrapolation relations, provided that 
the conserved variables at the points 3, 4, F, C, 7, 8, E, A, 10, H, 11, 12, R, S are known. The conserved 
variables at the cell centers labelled with letters are obviously known, while they must be interpolated at 
the points labelled with numbers. The point 3 is located between the cell centers H and I, the point 7 
between G and H , the point 10 between E and F. A more uncertain situation holds for the remaining 
points 4, 8, 11 and 12. Indeed, depending on the coarsening, points 4, 8 and 11 may be located between 
cell centers M and L (present situation) or between cell centers L and N, while point 12 may be located 
between Q and P or between P and O. A straightforward linear interpolation based on the values at the 
surrounding cell centers provides the values of the conserved variables at the points 3, 7 and 10. 

The uncertainty in the interpolation of these variables at point 4 can be overcome by using the 
following relation: 


<74 = <7l + 


VL ~ VM VN ~ VL 


(17) 
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where q indicates the conserved variables and y the ordinates, while the subscript letters correspond to 
the cell centers in Fig. 13. Furthermore, U4 = 1 if yp < yp, otherwise U4 = 0. Using Eq. (17), the 
conserved variables at point 4 are linearly interpolated using the values pertaining to the cell centers M 
and L if point 4 is located between these cell centers. Otherwise, the interpolation uses values at cell 
centers L and N. Similarly, at points 8, 11, 12, the following relations can be used: 
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Eq. (20), U 12 = 1 if yp < yp, otherwise U 12 = 0. 

Using the values at the appropriate cell centers and the interpolated values at the other points, the 
left and right states at the vertical cell edges 1, 5, 9, 13, 14 can be computed and the fluxes at these cell 
edges evaluated. 

At this point, we see from Fig. 13 that the topology at the right side of cell H is general, while this 
is not true for the topology on the left side. Indeed, the cell center E can also belong to a coarsening 
element with smaller cells located at the left side of E as illustrated in Fig. 14, where even the cell F is 
in such a situation. This would bring in additional topological possibilities and complicate the evaluation 
of the fluxes. This in turn is not necessary if the coarsening element evaluation is performed sweeping 
form the highest i value (corresponding to the outer boundary) to the smallest one (corresponding to 
the left boundary of the coarsening region downstream of the body). In such a way, the possible errors 
in the flux evaluation at the cell edges of the cell E would be corrected when computing the coarsening 
element including the cells A , B, E. Finally, it must be emphasized that the cell H has 5 cell edges with 
different fluxes. This situation must be taken into account when computing the residuals at this cell. 
The previous procedure with minor modifications may be utilized for the evaluation of the fluxes at the 
coarsening element cell edges for the other three coarsening regions. 


2.2 Far-Field Coarsening Results 

We considered three test cases corresponding to subsonic and transonic flow conditions about a NACA 
0012 airfoil. The first one is the transonic flow with an upstream Mach number = 0.80 at an angle of 
attack of 1.25°. Figure 15 presents the Mach number distribution on the airfoil surface (continuous lines) 
computed by applying the grid coarsening approach. In this figure, this distribution is compared with 
the original grid Mach number distribution (symbols). The two sets of results in Fig. 15 coincide. The 
comparison of the aerodynamic parameters enforces this finding. The coarsened grid computation gives a 
lift coefficient Cl = 0.3529 and a drag coefficient Cp = 0.0232, practically coincident with the original 
grid computation results: Cl = 0.3535, Cp = 0.0232. In this test case, the far-field coarsening reduces 
the original 17272 cells to only 7601 cells. Only 44% of the original cells are computed resulting in a 34% 
reduction in CPU time. The cell center distribution is represented in Fig. 16 for the original grid, while 
the effect of the far-field coarsening is evident in the cell center distribution in Fig. 17. The clustering 
of the grid near the body is also evident in the far-field in Fig. 16, while it is completely eliminated in 
Fig. 17, where the far-field mesh presents an aspect ratio not far from one. For comparison, this test case 
has been recomputed using an uncoarsened grid of 8366 cells with 7747 cells, located in the flow field, 
to be computed, i.e., a cell number slightly higher than the number of cells considered in the previous 
coarsened grid computation. The computed aerodynamic parameters are Cl = 0.3610, Cp = 0.0241, 
which are significantly different from the reference results: Cl = 0.3535 and Cp = 0.0232. 

The second test case is the subsonic flow with an upstream Mach number = 0.63 at an angle 
of attack of 2°. Figure 18 shows the Mach number distribution on the airfoil surface computed by 
applying the grid coarsening approach (continuous lines), compared with the original grid Mach number 
distribution (symbols). Fig. 18 again shows the agreement of the two distributions. In this test case, 
the original grid lift and drag coefficients are Cl = 0.3345, Cp = 0(10 -4 ) while the coarsened grid 
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computation gives Cl = 0.3341 and Cd = 0(10 -4 ). Furthermore, only 40% of the original cells are 
computed and the CPU time is reduced to 47%. 

Finally, Fig. 19 shows the Mach number distributions for the last test case, transonic flow with 
an upstream Mach number = 0.85 at an angle of attack of l.°, and once again, indicates their 
agreement. The computed original grid aerodynamic parameters are Cl = 0.3726, Cd = 0.0588 while 
the coarsened grid results are Cl = 0.3721 and Cd = 0.0589. For this test case, only 50% of the 
original cells are computed and the CPU time is reduced to 38%. All the computed test cases prove that 
the far-held coarsening practically does not alter the computed results, while significantly reducing the 
computational time. One can easily imagine even more dramatic results for three-dimensional cases. 

3. INVERSE DESIGN OF AIRFOILS 


3.1 Objective Function 

The inverse design of airfoils consists of finding the body geometric shape whose pressure distribution 
along the body surface, p{s,£), matches a target pressure distribution, p{s), where s is the curvilinear 
coordinate measured along the surface. A discrete objective function, I, may be defined as: 

1 Nb 

I (0 = yw 5 >' /'']"• ( 21 ) 

i = 1 


where the body length has been subdivided into Nb intervals, whose center locations are s,: (i = 1, . . . , Nb). 
The discrete interpolated pressures are p, = p{si,£) and the discrete target pressures are pi = p(s,:), where 
£ represents the design parameters. 

3.2 Adjoint Formulation 

An auxiliary flow solver is used only for the computation of smooth sensitivity derivatives. It should 
be characterized by added artificial viscosity in order to smooth the shock. The simplest approach which 
will meet our goals is to use a central scheme with an adequate level of artificial viscosity, as described 
in Ref. 21. This auxiliary flow solver is used only in the discrete adjoint formulation described below. 

The technique presented in Ref. 21 brings us to the following equation for the computation of the 
derivatives of the objective function: 


dl 

dij 



( 22 ) 


where the adjoint variables A*, are evaluated according to the following conditions: 


<9A 

dt 




(23) 


for i = 1, . . . , N x M and j = 1, . . . , n, with n the number of design parameters, and repeated indices 
summed over their range. In the previous equations, u, corresponds to the vector of the physical variables 
at each cell center, i: 


u, : = [ p , p, u, v]f, (24) 

while w*(u,£) = 0 represents the discretized system of governing partial differential equations corre- 
sponding to the auxiliary flow solver. 

Equation (23) represents a set of 4 x N x M adjoint equations. In addition, the terms (dw^/du , : )^ 
at the cell centers away from the body can be evaluated directly from the system of discretized flow 
equations, w* = 0, and depend on the solution of the flow field. At the cell centers close to the body, they 
are numerically evaluated following the procedure suggested in Ref. 25, which determines these values 
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by perturbing m at a cell center close to the body and by determining the corresponding variations 
of the residuals at the considered cell center and the neighboring cell centers. Then (dw^/duj)^ is 
approximated by the ratio between the variations of the residuals and the corresponding variations of 
In this computation, the corresponding changes of the boundary conditions are also taken into account. 
The terms (<9w k/d£j) u are numerically evaluated in a similar way by perturbing the design parameters 

by evaluating the corresponding variations of the fluid dynamic properties at the row of cell centers 
surrounding the body on the interior and by determining the corresponding changes of the residuals at 
the cell centers surrounding the body and located in the flow field. The terms (dl/d u,)^ can be easily 
evaluated from the objective function expression, while the terms (dl/d£j) u are evaluated by perturbing 
the design parameters and by evaluating the corresponding changes of the objective function. 

The fluid dynamic variables computed with the accurate, basic flow solver are used to compute the 
terms (dl/d u,;)^ and (dw^/duj)^. Then, the adjoint equations are solved for the 4 x N x M adjoint 
variables, Ak- During this computation, the far field boundaries are moved closer to the body with 
respect to their location for the flow field computation. This approximation allows to significantly reduce 
the CPU time required to evaluate the adjoint variables without significantly change their values at the 
cell centers located in the vicinity of the body. During the adjoint equation solution, the adjoint variables 
corresponding to the cell centers located inside the body are set to zero and never changed during the 
computation. Finally, Eq. (22) is used to evaluate the sensitivity derivatives of the objective function, 

M/dti- 

At this point we faced the choice of using the coarsened grid or resorting to an optimization process 
on the complete grid. We decided to adopt the second choice and to use the complete uncoarsened grid 
because of the following reasons: i) the domain employed to solve the adjoint equations is significantly 
reduced with respect to the one employed to solve the flow equations, so that the savings in cell center 
evaluation would have been very small; ii) the coarsening requires a special treatment at the de-refinement 
elements. The computational cost of this treatment may not be small, in comparison to the CPU saving 
due to the reduced number of cell centers to be computed, if few iblanked cell centers are involved. We 
concluded that the possible saving would have been very small, if any. Thus, before transferring the 
fluid dynamic data to the optimizer, we had to evaluate these data at all cell centers corresponding to 
the original uncoarsened grid. This problem was simply solved by interpolating the solution from the 
coarsened grid to the original one by a linear interpolation based on the values computed at the active 
cell centers. 


3.3 Progressive Optimization 

The goal of the progressive convergence technique is to determine the minimum of the objective 
function, while progressively reaching a converged solution for the flow field and the adjoint equations. 
The current implementation of the progressive algorithm is as follows: 

1. Define the coarsest mesh and several finer levels. Each finer mesh is obtained from the corresponding 
coarser mesh by doubling the number of intervals in each direction. 

2. Start with an initial set of design variables. 

3. Start the flow computations on the coarsest grid. 

4. Advance the flow solver for several iterations. 

5. Evaluate the adjoint equations by advancing the adjoint solver for several iterations, with the 
condition that the residuals drop at least one order of magnitude with respect to the value it 
assumes at the beginning of the procedure. 

6. Compute the objective function gradient, V/. 

7. Update the design variables according to the relation: 


( 1 + 1 = c' _ a — 
^ ' d£j' 


( 25 ) 


where aj are positive parameters. 
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8. Repeat steps 4 to 7 until the gradient of the objective function is sufficiently decreased (from half 
to one order of magnitude), with the condition that the residuals of the flow and of the adjoint 
equations drop progressively. 

9. Move to the corresponding finer mesh and interpolate the computed flow solution to the finer grid. 

10. Repeat steps 4 to 7 until the gradient of the objective function is decreased half to one order of 
magnitude more. 

11. Repeat steps 9 and 10 until the finest grid is reached. 

12. Repeat steps 4 to 7 until the objective function has sufficiently decreased. 


The coefficient a.j in Eq. (25) is evaluated as a,j = b cj. The coefficient b represents the minimum 
value of cij and is given by: 


b 


\Mc\ 

k \VjI\Ur 


(26) 


where A£ c is a typical (order of magnitude) change of the design parameters, A: is a constant ranging 
from 40 to 100 and depending on the test case, and \YjI\^ nax is the largest absolute value assumed by 
the sensitivity derivatives of the objective function during the first global iterations through steps 2 to 
6. The term Cj is an amplification factor with a minimum value equal to 1 and a maximum value cm 
which depends on the convergence of the normalized gradient of the objective function. In particular, cm 
is equal to 10 if the convergence is lower than one order of magnitude and is equal to 40 if it is higher 
than two orders, while it varies linearly between 10 and 40 for intermediate values of the convergence of 
the normalized gradient of the objective function. At the beginning of the computations, cj is set to 1, 
then it is increased by 50% if the corresponding sensitivity derivative of the objective function maintains 
its sign, while it is decreased by 50% if the sign changes. This approach allows large changes for the 
design variables whose sensitivity derivatives maintain the same sign, while it assigns small changes to 
design variables whose sensitivity derivatives are changing their sign. We have devised this technique 
because of the very slow convergence of the method of steepest descent, corresponding to Cj equal to 
one. Our approach of using a smoothed adjoint solver based on an auxiliary flow solver leads to an 
approximate gradient evaluation, which caused convergence problems with classical optimizers. The 
use of the amplification factor cj described above, significantly improves convergence, compared to the 
method of steepest descent. 

In the computed results, we have used three grid levels and we have considered the objective function 
sufficiently decreased when the norm of its gradient has decreased 3 orders of magnitude with respect to 
the value it assumes after the first global iteration through steps 2 to 6. 


3.4 Inverse-Design Results 

The inverse design of airfoils in subsonic and transonic flow conditions has been considered. Following 
the procedure used in Ref. 26, the airfoil surface shape is described as: 

N a 

Y(x/c) = Y^ZiUx/c), (27) 

i=i 


where Y = y/c , with y being the airfoil ordinate and c the airfoil chord length. The shape functions 
))(.r/c) are the ordinates of the N a specified airfoils, which depend on the non-dimensional abscissa x/c. 
Finally, the weighting terms are the design parameters. 

First, we considered the subsonic flow about an NACA 0012 airfoil with an undisturbed Mach number 
Moo = 0.63 at an angle of attack of 2.°, with the aim to check the ability of the method to bring the 
convergence of the optimization process to machine zero. The reference flow solution on the finest mesh 
converged to machine zero was used to define the pressure distribution on the airfoil surface. The airfoil 
optimization test was performed with only one design variable. The target airfoil is the NACA 0012 
airfoil so that the target value of the design parameter is £ = 1. We assumed an initial airfoil shape 
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corresponding to a value of the design parameter £ = 2.5. The optimization process was performed using 
three mesh levels: 32 x 29, 64 x 59, 128 x 119 cells, with the larger number of cells being located in the x 
direction. A line of cell centers was located at y = 0 at all mesh levels. Fig. 20 shows the convergence 
history of the logarithm of the objective function, log 10 (7), and of the logarithm of the magnitude of the 
normalized gradient of the objective function, log 10 |V7|. The abscissa in Fig. 20 is the required work, 
with a unit of work considered as the computational time required to run a single analysis solution to 
convergence on the finest mesh. Fig. 20 shows that log 10 (7) converged to 10 -22 and log 10 |V7| converged 
to 10 -14 . Meanwhile, the design parameter reached its target value, equal to 1, with a residual error of 
10 _u , These results indicate that the optimization process reached the machine zero. 

Then, we considered the transonic flow about an NACA 0012 airfoil with an undisturbed Mach number 
Mqo = 0.80 at an angle of attack of 1.25°. The airfoil optimization test was performed with N a = 4 and 
lj-l ’4 as the shape of four specified airfoils: NACA 0012, NACA 64i-412, NACA 64 2 A215, NACA 65o- 
415. The target airfoil is the NACA 0012 airfoil so that the target set of coefficients is £ = (1, 0, 0, 0). 
We assumed an initial airfoil shape corresponding to the set of coefficients £ = (1.5, 0.5, 0.5, 0.5). The 
optimization process was performed using the three mesh levels employed in the previous test case. Fig. 21 
shows the convergence history of the logarithm of the objective function, log 10 (7), and of the logarithm of 
the magnitude of the normalized gradient of the objective function, log 10 |V7|. The results indicate that 
the work required to converge the gradient of the objective function by 3 orders of magnitude is equal to 
1.07. In our progressive optimization strategy, large gains in efficiency are obtained because most of the 
geometrical changes occur at the coarsest mesh level at a very low computational cost. Fig. 22 shows the 
target pressure distribution on the airfoil surface (continuous line) and the converged pressure distribution 
(symbols). The two sets of results are practically coincident even in the very sensitive shock transition 
region, proving that the optimized solution is sufficiently converged. Fig. 23 presents the target airfoil 
shape (continuous line), the initial airfoil shape (dashed line), and the optimized airfoil shape (symbols). 
The figure shows that the optimized airfoil shape practically coincides with its target value although the 
optimization process starts from an airfoil shape very different from the optimized one. 

Finally, we repeated the previous optimization test case by changing only the flow conditions. Here 
we considered the sometime difficult transonic flow about an NACA 0012 with an undisturbed Mach 
number = 0.85 at an angle of attack of l.°. Fig. 24 shows the convergence history of the logarithm 
of the objective function, log 10 (7), and of the logarithm of the magnitude of the normalized gradient of 
the objective function, log 10 |V7|. The results indicate that the work required to converge the gradient 
of the objective function by 3 orders of magnitude is equal to 2.33. Fig. 25 shows the target pressure 
distribution on the airfoil surface (continuous line) and the converged pressure distribution (symbols). 
The two sets of results are practically coincident even in the very sensitive shock transition region, proving 
that the optimized solution is sufficiently converged. 

CONCLUSIONS 


The first part of this report reviews a new approach to the computation of inviscid flows using Carte- 
sian grid methods. The crux of the method is the curvature-corrected symmetry technique (CCST) 
developed by the present authors for body-fitted grids. The method introduces ghost cells near the 
boundaries whose values are developed from an assumed flow-field model in vicinity of the wall, consist- 
ing of a vortex flow with a symmetric distribution of entropy and total enthalpy along a local normal to 
the body surface. This model satisfies the normal momentum equation, the non-penetration boundary 
condition and Crocco’s relation. The CCST boundary condition was shown to be substantially more 
accurate than traditional boundary condition approaches. This improved boundary condition has been 
adapted to a Cartesian mesh formulation, the Ghost Body-Cell Method. In this method, all cell centers 
exterior to the body are computed with fluxes at the four surrounding cell edges. There is no need for spe- 
cial treatment corresponding to cut cells which complicate other Cartesian mesh methods. The approach 
does require surface normals and radii of curvature. However, since the body is defined pointwise and is 
independent of the mesh, it may be easily defined by a sufficient point density to allow straightforward 
accurate numerical calculation of the surface normals and radii of curvature. 

The investigation of the compressible inviscid flow about a circular cylinder has shown that the non- 
penetration condition is satisfied in the limit of vanishing cell size, with errors in the velocity component 
normal to the body of the order of 10 -4 for the finest grid employed (with free-stream velocities scaled 
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to 1). The Cartesian grid results have then been compared with the results computed using the classical 
second-order pressure extrapolation technique and the curvature-corrected symmetry technique on body 
fitted polar grids. The computed contours of total pressure have been compared. This result proves 
that GBCM on Cartesian grids presents dramatic advantages with respect to the widely used second- 
order pressure extrapolation technique on body-fitted grids. A mesh refinement study has indicated that 
GBCM is second-order accurate in space. 

Results for the transonic flow about an NACA 0012 airfoil have shown that GBCM can be easily 
extended to practical body shapes. The computed Mach number distributions compare quite well with 
body-fitted grid results obtained using the CCST technique. The computed lift and drag coefficients 
match published values considered to be accurate. Importantly, the results have been conclusively shown 
to be independent of the position of the leading edge with respect to the grid. The flexibility of the 
method has been finally illustrated with an application to the flow over a multiple-body configuration, 
the bi-NACA configuration. Comparison with published results has further confirmed the accuracy of 
the present results. 

The original version of the GBCM, as with all Cartesian mesh formulations, has the property that any 
grid clustering near the body is maintained to the far field boundary. This is not a significant penalty in 
2D calculations, but it may become significant in 3D calculations. In the second part of this report, we 
have introduced a far-field grid coarsening to address this issue. The technique we have adopted is based 
on iblanking of cell centers. The computational domain has been subdivided in five regions, the central 
one including the body and characterized by no coarsening and four external regions, located above, 
below, upstream and downstream of the body where the grid coarsening is applied. Only the coalescence 
of two smaller cells into a larger one is allowed, so that the coarsening is essentially anisotropic. The 
coarsening is applied only in the direction which is characterized by a smaller grid spacing and tends to 
determine near-square cells. The special treatment required at the connection between two smaller cells 
and a larger one determined by the coalescence of two smaller original cells is described in detail. 

The far-field grid coarsening has been applied to subsonic and transonic flows about an NACA 0012 
airfoil, indicating that the results are practically unaffected by the coarsening. Meanwhile, the number 

of the computed cells is reduced to 40% 50% and the CPU time to 34% 47% of the original CPU 

time, depending on the particular test case. 

The Ghost Body-Cell Method on Cartesian far-field coarsened grids appears to be particularly well- 
suited for design optimization problems. The third part of the report addresses this issue and presents 
an efficient formulation for the robust inverse design of compressible fluid flow problems. The design 
approach has five essential features: First, the procedure utilizes the accurate GBCM flow solver on 
coarsened Cartesian grids. Second, approximate but efficient design sensitivities are obtained using a 
discrete adjoint formulation. The adjoint problem is based on an artificially dissipative flow solver in 
order to obtain robust sensitivity derivatives in the presence of noisy or non-smooth objective functions. 
This adjoint solver is only used for the solution of robust sensitivities. The next feature of the procedure is 
the use of quasi-time-dependent adjoint equations. The fourth feature involves what we term progressive 
optimization , whereby a sequence of operations, containing a partially converged flow solution, followed 
by a partially converged adjoint solution followed by an optimization step is performed. The last feature 
is the use of progressively finer grids for the progressive solution of the flow field, while only using a 
reduced domain for the solution of the adjoint equations. 

The inverse design of airfoils in transonic flow conditions indicated that the methodology is accurate, 
robust and highly efficient, with converged design optimizations produced in no more than the amount of 
computational work to perform two flow analysis on the finest mesh. In addition, the inverse design of a 
NACA 0012 airfoil in subsonic flow conditions showed that the procedure is able to bring the convergence 
of the optimization process to machine zero. 

Overall, we observe that the work to perform inverse design of airfoils on Cartesian grids using progres- 
sive optimization procedures is comparable to the work required by the design optimization procedure on 
body fitted grids, previously developed. Furthermore, the Cartesian grid design procedure only requires 
the generation of a simple Cartesian mesh which remains unchanged during the optimization iterations, 
while the body-fitted grid design procedure requires an ad hoc grid generation process and a continuous 
remeshing during the optimization process. This difference represents a significant advantage of the de- 
sign procedure on Cartesian grids and makes this procedure much better suited for complex geometry 
optimization. 
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Table 1: Effects of the leading edge position. 
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Figure 2: Cartesian grid: cell center net, interior cell centers close to the body (open symbols) and 
symmetrically reflected exterior cell centers (solid symbols). 
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Figure 3: Double-valued ghost cell centers near a sharp trailing edge. 
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Figure 4: Grid convergence of flow over circular cylinder using Ghost Body-Cell Method. Lo-norm of the 
surface normal velocity component and Mach number at the leading and trailing edges. 
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Figure 5: Isobars of total pressure. Ghost Body-Cell Method on Cartesian grid. 
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Figure 7: Surface Mach number distribution for NACA 0012 airfoil at = 0.80, a = 1.25°. 

Cartesian grid GBCM (solid line) and body-fitted CCST (symbols). 
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Figure 8: Cartesian grid GBCM calculation of NACA 0012 airfoil at = 0.85, a = 0°. Surface 
Mach number distribution in the leading edge region. Calculations performed with various leading-edge 
locations with respect to the grid as indicated by the inset. 
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Figure 15: Surface Mach number distribution for NACA 0012 airfoil at = 0.80, a 

Coarsened grid (solid line) and uncoarsened grid (symbols). 
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Figure 16: ,\l x = 0.80, a 


1.25°. Cell center net for the original uncoarsened grid. 
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Figure 17: = 0.80, a 


1.25°. Cell center net for the grid with far-held coarsening. 
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Figure 18: Surface Mach number distribution for NACA 0012 airfoil at = 0.63, a = 
grid (solid line) and uncoarsened grid (symbols). 
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Figure 19: Surface Mach number distribution for NACA 0012 airfoil at = 0.85, a = 1°. Coarsened 
grid (solid line) and uncoarsened grid (symbols). 
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Figure 20: Convergence history for inverse design of a subsonic NACA 0012 airfoil. = 0.63, a = 2° 
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Figure 21: Convergence history for inverse design of a transonic NACA 0012 airfoil. = 0.80, 

a = 1.25°. 
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Figure 22: Transonic NACA 0012 airfoil. = 0.80, a = 1.25°. Target (continuous line) and 

optimized (symbols) pressure distribution on airfoil surface. 
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Figure 24: Convergence history for the inverse design of a transonic NAC'A 0012 airfoil. = 0.85, 

OL — 1 ° 
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